Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

lprior argument to prior, enabling separate lprior variables #1724

Open
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

n-kall
Copy link

@n-kall n-kall commented Jan 14, 2025

This would allow for tagging priors for inclusion in different lprior variables, enabling priorsense to selectively power-scale priors. It would fix #1585

@n-kall n-kall changed the title initial work on lprior tags lprior argument to prior, enabling separate lprior variables Jan 14, 2025
@paul-buerkner
Copy link
Owner

paul-buerkner commented Jan 15, 2025

I think this looks pretty good already.

You mentioned that so far it only works with the b coefficents but I couldn't readily see this limitation in your code. Can you elaborate?

Also, could you show an example model here to illustrate the features use?

Thank you for working on this issue!

@paul-buerkner paul-buerkner added this to the brms 2.23.0 milestone Jan 15, 2025
@n-kall
Copy link
Author

n-kall commented Jan 15, 2025

Certainly, here is an example:

Basically, the current code works for priors that are set on the coefficient level, i.e. has_coef_prior == TRUE see here

I wasn't sure how to properly handle the lprior tagging for a class when there is no individual coef prior. Now the lprior variables are created, but not added to if there is no individual coefficient prior.

prior <- c(
  prior(normal(0, 5), coef = "K1", lprior = "covariates"), # works
  prior(normal(0, 1), coef = "N1", lprior = "treatment"), # works
  prior(normal(0, 100), class = "Intercept", lprior = "Intercept"), # does not work
  prior(normal(0, 10), class = "sigma", lprior = "sigma") # does not work
)


x <- brm(yield ~ N + K + P, data = npk,
         prior = prior, backend = "cmdstanr",
         empty = TRUE
         )

prior_summary(x)

s <- stancode(x)

s
##  > prior_summary(x)
##           prior     class coef group resp dpar nlpar lb ub     lprior       source
##          (flat)         b                                                  default
##    normal(0, 5)         b   K1                             covariates         user
##    normal(0, 1)         b   N1                              treatment         user
##          (flat)         b   P1                                        (vectorized)
##  normal(0, 100) Intercept                                   Intercept         user
##   normal(0, 10)     sigma                             0         sigma         user
// generated with brms 2.22.8
functions {
}
data {
   int<lower=1> N;  // total number of observations
   vector[N] Y;  // response variable
   int<lower=1> K;  // number of population-level effects
   matrix[N, K] X;  // population-level design matrix
   int<lower=1> Kc;  // number of population-level effects after centering
   int prior_only;  // should the likelihood be ignored?
 }
 transformed data {
   matrix[N, Kc] Xc;  // centered version of X without an intercept
   vector[Kc] means_X;  // column means of X before centering
   for (i in 2:K) {
     means_X[i - 1] = mean(X[, i]);
     Xc[, i - 1] = X[, i] - means_X[i - 1];
   }
 }
 parameters {
   vector[Kc] b;  // regression coefficients
   real Intercept;  // temporary intercept for centered predictors
   real<lower=0> sigma;  // dispersion parameter
 }
 transformed parameters {
   real lprior = 0;  // prior contributions to the log posterior
   real lprior_covariates = 0;
   real lprior_treatment = 0;
   real lprior_Intercept = 0;
   real lprior_sigma = 0;
   lprior += normal_lpdf(b[1] | 0, 1);
   lprior_treatment += normal_lpdf(b[1] | 0, 1);
   lprior += normal_lpdf(b[2] | 0, 5);
   lprior_covariates += normal_lpdf(b[2] | 0, 5);
   lprior += normal_lpdf(Intercept | 0, 100);
   lprior += normal_lpdf(sigma | 0, 10)
     - 1 * normal_lccdf(0 | 0, 10);
 }
 model {
   // likelihood including constants
   if (!prior_only) {
     target += normal_id_glm_lpdf(Y | Xc, Intercept, b, sigma);
   }
   // priors including constants
   target += lprior;
 }
generated quantities {
   // actual population-level intercept
   real b_Intercept = Intercept - dot_product(means_X, b);
 }

@paul-buerkner
Copy link
Owner

I see. okay, I will take a look at make it work for all priors even without the coef specified

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

Interface for defining which prior log densities are stored for priorsense
2 participants